
*Welfare anaylsis based on matlab results

clear all

*Import the assembled dataset:
use "$root/Data/Produced/EV_welfare04.dta", clear


order price height weight var_cost typkw ev hyrid diesel car_size size_hh2 size_hh3 size_hh4 size_hh5 age_kw2 age_kw3 ev_urban2 ev_urban3 ev_dist_ev ev_range5 ev_year2 ev_year3 ev_pv ev_ho ev_wealth2 ev_wealth3 ev_wealth4 brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE2 car_cat_FE3 car_cat_FE4 car_cat_FE5 car_cat_FE6 car_cat_FE7 car_cat_FE8 env_friendly resid price_wealth2 price_wealth3 price_wealth4, after(choice)


*For belief error to be properly calculated var_costs need to be zero, where they are missing: 

replace var_cost=0 if missing(var_cost)

sort id car_option

*Set up the matrices and vectors needed for the computation:
mata

ndraws=150

b_f= (0.021183645,-2.593121925,0.154611204,-0.321820168,0.037136405,0.167413791,0.378321531,0.614503527,0.858999103,-0.005445116,-0.010654642,0.136554812,-0.083769393,-0.022510516,0.007285076,0.06246904,1.476449254,2.260850432,0.424973893,0.853075696,0.6579357,1.4180835,0.027265514,0.980066634,0.090419384,0.171157248,-0.142510175,-0.661616614,-0.177476937,0.595728981,0.214934965,0.378106038,-0.080343758,0.394597341,0.734410779,0.078025015,1.019999496,0.291640034,-0.620171333,-0.126470184,-0.114005609,-0.169462474,0.335901857,-0.232550483,0.070071269,0.05306741,0.003621919,0.006019563,0.024778481)'
/*
** Define the estimated vectors of random coefficients:
*/
beta_f=J(1,ndraws,b_f)
rseed(1234)
beta_p = rnormal(1,ndraws,-0.080965469,0.013364931)
beta_he = rnormal(1,ndraws,0.263640717,0.000661958)
beta_we = rnormal(1,ndraws,0.000192522,0.000045627)
beta_var = rnormal(1,ndraws,-0.040212482,0.001878464)

beta= (beta_p \ beta_he \ beta_we \ beta_var \ beta_f)

mata drop beta_he beta_we b_f


wealth_group=.
st_view(wealth_group,.,"wealth_quart2 wealth_quart3 wealth_quart4")


/*
Define the price coefficients:
*/
b_p=(0.003621919,0.006019563,0.024778481)'



b_w=wealth_group*b_p

mata drop wealth_group b_p



b_price=J(rows(b_w),1,beta_p) + J(1,ndraws,b_w)

gamma=beta_var:/b_price

mata drop beta_p b_w beta_var

end

timer on 1
mata
X=.
st_view(X,.,"price height weight var_cost typkw ev hyrid diesel car_size size_hh2 size_hh3 size_hh4 size_hh5 age_kw2 age_kw3 ev_urban2 ev_urban3 ev_dist_ev ev_range5 ev_year2 ev_year3 ev_pv ev_ho ev_wealth2 ev_wealth3 ev_wealth4 brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE2 car_cat_FE3 car_cat_FE4 car_cat_FE5 car_cat_FE6 car_cat_FE7 car_cat_FE8 env_friendly resid price_wealth2 price_wealth3 price_wealth4")


l_fit=exp(X*beta)

var_cost=.
st_view(var_cost,.,"var_cost")



/*
for (i=1; i<=ndraws; i++) {
st_addvar("double", "fit_" + strofreal(i))
st_store(.,"fit_" + strofreal(i),l_fit[,i])
}
*/
id=st_data(.,"id")
V=panelsetup(id, 1)

    ncol   = ndraws
    bigtot = J(rows(id), ncol, .)
    for (i = 1; i <= rows(V); i++) {
        X1 = panelsubmatrix(l_fit, i, V)
        bigtot[|V[i,1], 1 \ V[i, 2], ncol|] = J(rows(X1), 1, colsum(X1))
}

mata drop X1 ncol V

_editmissing(l_fit,0)

prob_all=l_fit:/bigtot
mata drop l_fit X 

prob=mean(prob_all')

be=prob_all :* (J(1, cols(gamma), var_cost)-(gamma:*var_cost))

be_m=mean(be')

gam=mean(gamma')

mata drop prob_all var_cost be

st_addvar("double", "prob")
st_store(.,"prob",prob')

st_addvar("double", "be_init")
st_store(.,"be_init",be_m')

st_addvar("double", "gamma")
st_store(.,"gamma",gam')

CS_init=mean((((-1):/b_price):*log(bigtot))')


mata drop bigtot

st_addvar("double", "CS_init")
st_store(.,"CS_init",CS_init')


mata drop prob CS_init be_m
end
timer off 1
timer list 1



**Set the directory to store the results:
cd "$root/Results/estimation_results/counterfactuals"


*First calculate the predicted probabilites for each category:

frame put prob et_co2 average_use et_verbrauch el_verbrauch car_tax pv_tax drivetype car_option prhnewprice, into(init_prob)

frame change init_prob

gcollapse (mean) prob et_co2 average_use et_verbrauch el_verbrauch car_tax pv_tax prhnewprice (firstnm) drivetype, by(car_option)


*We use the number of 9230 cars as the average registrations per year

gen co2_ann=prob*9230*et_co2*average_use/1000
gen emb_co2_ann=0
replace emb_co2_ann=prob*9230*average_use/100*et_verbrauch*0.523 if drivetype==1
replace emb_co2_ann=prob*9230*average_use/100*et_verbrauch*0.445 if drivetype==2
replace emb_co2_ann=prob*9230*average_use/100*el_verbrauch*0.139 if drivetype==3
replace emb_co2_ann=prob*9230*average_use/100*el_verbrauch*0.139 + prob*9230*average_use/100*et_verbrauch*0.523 if drivetype==4


gen petrol_usage=prob*9230*average_use/100*et_verbrauch
gen petrol_tax=petrol_usage*0.7312
replace petrol_tax=petrol_usage*0.7587 if drivetype==2

gen tax_total=prob*9230*car_tax
gen tax_pv=prob*9230*pv_tax


gen veh_tax=prob*9230*prhnewprice*0.04
replace veh_tax=0 if drivetype==3

gcollapse (sum) prob co2_ann emb_co2_ann petrol_tax tax_total tax_pv veh_tax, by(drivetype)

export excel "init_prob.xlsx", firstrow(variables) replace

frame change default
frame drop init_prob


*Also calculate the predicted probabilites for each category by iwealth group:
frame put prob et_co2 average_use et_verbrauch el_verbrauch car_tax pv_tax  drivetype car_option wealth_group prhnewprice, into(init_prob_inc)
frame change init_prob_inc

gcollapse (mean) prob et_co2 average_use et_verbrauch el_verbrauch car_tax pv_tax prhnewprice (firstnm) drivetype, by(car_option wealth_group)

gen co2_ann=prob*2307.5*et_co2*average_use/1000

gen emb_co2_ann=0
replace emb_co2_ann=prob*2307.5*average_use/100*et_verbrauch*0.523 if drivetype==1
replace emb_co2_ann=prob*2307.5*average_use/100*et_verbrauch*0.445 if drivetype==2
replace emb_co2_ann=prob*2307.5*average_use/100*el_verbrauch*0.139 if drivetype==3
replace emb_co2_ann=prob*2307.5*average_use/100*el_verbrauch*0.139 + prob*2307.5*average_use/100*et_verbrauch*0.523 if drivetype==4


gen petrol_usage=prob*2307.5*average_use/100*et_verbrauch
gen petrol_tax=petrol_usage*0.7312
replace petrol_tax=petrol_usage*0.7587 if drivetype==2


gen tax_total=prob*2307.5*car_tax
gen tax_pv=prob*2307.5*pv_tax

gen veh_tax=prob*2307.5*prhnewprice*0.04
replace veh_tax=0 if drivetype==3

gcollapse (sum) prob co2_ann emb_co2_ann petrol_tax tax_total tax_pv veh_tax, by(drivetype wealth_group)

export excel "init_prob_byinc.xlsx", firstrow(variables) replace

frame change default
frame drop init_prob_inc


******************************
**COUNTERFACTUAL:
**Introduction of a 4'000 CHF subsidy for EV analog to Canton of Ticino
**

gen price_alt_tot=prhnewprice
replace price_alt_tot=prhnewprice-4000 if drivetype==3

gen price_alt=price_alt_tot/1000

forvalues a=2/4 {
	
	gen price_w_alt`a'=price_alt * wealth_quart`a'

}

drop price_alt_tot

sort id car_option

mata
X=.
st_view(X,.,"price_alt height weight var_cost typkw ev hyrid diesel car_size size_hh2 size_hh3 size_hh4 size_hh5 age_kw2 age_kw3 ev_urban2 ev_urban3 ev_dist_ev ev_range5 ev_year2 ev_year3 ev_pv ev_ho ev_wealth2 ev_wealth3 ev_wealth4 brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE2 car_cat_FE3 car_cat_FE4 car_cat_FE5 car_cat_FE6 car_cat_FE7 car_cat_FE8 env_friendly resid price_w_alt2 price_w_alt3 price_w_alt4")


l_fit=exp(X*beta)

var_cost=.
st_view(var_cost,.,"var_cost")



/*
for (i=1; i<=ndraws; i++) {
st_addvar("double", "fit_" + strofreal(i))
st_store(.,"fit_" + strofreal(i),l_fit[,i])
}
*/
id=st_data(.,"id")
V=panelsetup(id, 1)

    ncol   = ndraws
    bigtot = J(rows(id), ncol, .)
    for (i = 1; i <= rows(V); i++) {
        X1 = panelsubmatrix(l_fit, i, V)
        bigtot[|V[i,1], 1 \ V[i, 2], ncol|] = J(rows(X1), 1, colsum(X1))
}

mata drop X1 ncol V

_editmissing(l_fit,0)

prob_all=l_fit:/bigtot
mata drop l_fit X 

prob=mean(prob_all')

be=prob_all :* (J(1, cols(gamma), var_cost)-(gamma:*var_cost))

be_m=mean(be')

mata drop prob_all var_cost be


st_addvar("double", "prob_count2")
st_store(.,"prob_count2",prob')

st_addvar("double", "be_count2")
st_store(.,"be_count2",be_m')

CS=mean((((-1):/b_price):*log(bigtot))')


mata drop bigtot

st_addvar("double", "CS_count2")
st_store(.,"CS_count2",CS')


mata drop prob CS be_m

end



*First calculate the predicted probabilites for each category:
frame put prob_count2 et_co2 average_use et_verbrauch el_verbrauch car_tax pv_tax drivetype car_option prhnewprice, into(count2)

frame change count2

gcollapse (mean) prob_count2 et_co2 average_use et_verbrauch el_verbrauch car_tax pv_tax prhnewprice (firstnm) drivetype, by(car_option)


*We use the number of 9230 cars as the average registrations per year

gen co2_ann=prob_count2*9230*et_co2*average_use/1000
gen emb_co2_ann=0
replace emb_co2_ann=prob_count2*9230*average_use/100*et_verbrauch*0.523 if drivetype==1
replace emb_co2_ann=prob_count2*9230*average_use/100*et_verbrauch*0.445 if drivetype==2
replace emb_co2_ann=prob_count2*9230*average_use/100*el_verbrauch*0.139 if drivetype==3
replace emb_co2_ann=prob_count2*9230*average_use/100*el_verbrauch*0.139 + prob_count2*9230*average_use/100*et_verbrauch*0.523 if drivetype==4


gen petrol_usage=prob_count2*9230*average_use/100*et_verbrauch
gen petrol_tax=petrol_usage*0.7312
replace petrol_tax=petrol_usage*0.7587 if drivetype==2


gen tax_total=prob_count2*9230*car_tax
gen tax_pv=prob_count2*9230*pv_tax
gen paid_subsidy=prob_count2*9230*4000 if drivetype==3
gen veh_tax=prob_count2*9230*prhnewprice*0.04
replace veh_tax=0 if drivetype==3


gcollapse (sum) prob_count2 co2_ann emb_co2_ann petrol_tax tax_total tax_pv paid_subsidy veh_tax, by(drivetype)

export excel "alt2_prob.xlsx", firstrow(variables) replace

frame change default
frame drop count2

*Also calculate the predicted probabilites for each category by income quartile
frame put prob_count2 et_co2 average_use et_verbrauch el_verbrauch car_tax pv_tax drivetype car_option wealth_group prhnewprice, into(count2_inc)
frame change count2_inc

gcollapse (mean) prob_count2 et_co2 average_use et_verbrauch el_verbrauch car_tax pv_tax prhnewprice (firstnm) drivetype, by(car_option wealth_group)

gen co2_ann=prob_count2*2307.5*et_co2*average_use/1000
gen paid_subsidy=prob_count2*2307.5*4000 if drivetype==3
gen tax_total=prob_count2*2307.5*car_tax
gen tax_pv=prob_count2*2307.5*pv_tax

gen emb_co2_ann=0
replace emb_co2_ann=prob_count2*2307.5*average_use/100*et_verbrauch*0.523 if drivetype==1
replace emb_co2_ann=prob_count2*2307.5*average_use/100*et_verbrauch*0.445 if drivetype==2
replace emb_co2_ann=prob_count2*2307.5*average_use/100*el_verbrauch*0.139 if drivetype==3
replace emb_co2_ann=prob_count2*2307.5*average_use/100*el_verbrauch*0.139 + prob_count2*2307.5*average_use/100*et_verbrauch*0.523 if drivetype==4


gen petrol_usage=prob_count2*2307.5*average_use/100*et_verbrauch
gen petrol_tax=petrol_usage*0.7312
replace petrol_tax=petrol_usage*0.7587 if drivetype==2
gen veh_tax=prob_count2*2307.5*prhnewprice*0.04
replace veh_tax=0 if drivetype==3


gcollapse (sum) prob_count2 co2_ann emb_co2_ann petrol_tax paid_subsidy tax_total tax_pv veh_tax, by(drivetype wealth_group)

export excel "alt2_prob_byinc.xlsx", firstrow(variables) replace


frame change default
frame drop count2_inc


**Calculate the change in consumer welfare 

*Decision consumer surplus change:
gen CS_change2=CS_count2-CS_init

*Change in belief error: 
gen be_change2=be_count2-be_init

frame put CS_change2 CS_init be_change be_init wealth_group id, into(cs2)
frame change cs2

gcollapse (firstnm) CS_change2 CS_init wealth_group (sum) be_change2 be_init  , by(id)
gcollapse (sum) CS_change2 CS_init be_change2 , by(wealth_group)
gen CS_perc_change2=CS_change2/abs(CS_init)

replace CS_change2=2307.5/5768.5*CS_change2

replace be_change2=2307.5/5768.5*be_change2

export excel "CS2.xlsx", firstrow(variables) replace

frame change default
frame drop cs2


****
*Evaluate the goodness of fit of the predicted probabilities: 
***

frame put id car_option prob choice drivetype wealth_group, into(prediction)
frame change prediction

gegen max_prob=max(prob), by(id)
gen predict = (max_prob==prob)

preserve
gcollapse (mean) prob (firstnm) drivetype (sum) choice, by(car_option wealth_group)


gcollapse (sum)prob choice, by(drivetype wealth_group)

gen predicted=round(5768.5*prob)

drop prob
rename choice observed

reshape wide predicted observed, i(wealth_group) j(drivetype)

export excel "average_predictions.xlsx", firstrow(variables) replace


restore

*Check stats overall for each category seperately: 
preserve
gcollapse (mean) prob (firstnm) drivetype (sum) choice, by(car_option)


gcollapse (sum)prob choice, by(drivetype)

gen predicted=round(23074*prob)

drop prob
rename choice observed


export excel "average_predictions_overall.xlsx", firstrow(variables) replace


restore


frame change default
frame drop prediction 

** -------------------------------------------------------------
** Repeat the same exercise for different level of subsidies
** 


*At level zero, the results represent the initial probability analysis: 
frame create init
frame change init

import excel "$root/Results/estimation_results/counterfactual/init_prob_byinc.xlsx", sheet("Sheet1") firstrow clear

rename prob prob_subs
gen subsidy=0
gen paid_subsidy=0


frame put _all, into(frame1)

frame change default
frame drop init

*At level zero the CS change is also zero
frame create cons1
frame change cons1

set obs 1 
gen CS_subs_pchange=0
gen subsidy=0
gen CS_subs_change=0
gen be_subs_change=0
frame change default	


drop price_alt price_w_alt2 price_w_alt3 price_w_alt4

local count=2


forvalues x=300(300)9900 {

timer on 1
	
	
	gen subsidy=`x'
	
	gen price_alt_tot=prhnewprice
replace price_alt_tot=prhnewprice-subsidy if drivetype==3

gen price_alt=price_alt_tot/1000

forvalues a=2/4 {
	
	gen price_w_alt`a'=price_alt * wealth_quart`a'

}



drop price_alt_tot

sort id car_option

mata {
X=. ;
st_view(X,.,"price_alt height weight var_cost typkw ev hyrid diesel car_size size_hh2 size_hh3 size_hh4 size_hh5 age_kw2 age_kw3 ev_urban2 ev_urban3 ev_dist_ev ev_range5 ev_year2 ev_year3 ev_pv ev_ho ev_wealth2 ev_wealth3 ev_wealth4 brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE2 car_cat_FE3 car_cat_FE4 car_cat_FE5 car_cat_FE6 car_cat_FE7 car_cat_FE8 env_friendly resid price_w_alt2 price_w_alt3 price_w_alt4") ;


l_fit=exp(X*beta) ;

var_cost=. ;
st_view(var_cost,.,"var_cost") ;

/*
for (i=1; i<=ndraws; i++) {
st_addvar("double", "fit_" + strofreal(i))
st_store(.,"fit_" + strofreal(i),l_fit[,i])
}
*/
id=st_data(.,"id") ;
V=panelsetup(id, 1) ;

    ncol   = ndraws ;
    bigtot = J(rows(id), ncol, .) ;
    for (i = 1; i <= rows(V); i++) { ;
        X1 = panelsubmatrix(l_fit, i, V) ;
        bigtot[|V[i,1], 1 \ V[i, 2], ncol|] = J(rows(X1), 1, colsum(X1)) ;
} ;

/*
mata drop X1 ncol V
*/
X1=. ;
V=. ;
_editmissing(l_fit,0) ;

prob_all=l_fit:/bigtot ;
/*
mata drop l_fit X 
*/
l_fit=. ;
X=. ;

prob=mean(prob_all') ;

be=prob_all :* (J(1, cols(gamma), var_cost)-(gamma:*var_cost)) ;

be_m=mean(be') ;

/*
mata drop prob_all
*/
prob_all=. ;
be=. ;

st_addvar("double", "prob_subs") ;
st_store(.,"prob_subs",prob') ;

st_addvar("double", "be_subs") ;
st_store(.,"be_subs",be_m') ;

CS=mean((((-1):/b_price):*log(bigtot))');

/*
mata drop bigtot
*/
bigtot=.;
be_m=. ;
st_addvar("double", "CS_subs");
st_store(.,"CS_subs",CS');

/*
mata drop prob CS
*/
prob=.;
CS=.;
}


*Create a frame to calculate the consumer welfare change:
frame put CS_subs CS_init be_subs be_init id wealth_group subsidy, into(cons`count')

frame put prob_subs et_co2 average_use et_verbrauch el_verbrauch car_tax pv_tax drivetype subsidy prhnewprice car_option wealth_group, into(frame`count') 
frame change frame`count'

gcollapse (mean) prob_subs et_co2 average_use et_verbrauch el_verbrauch car_tax pv_tax prhnewprice (firstnm) drivetype subsidy, by(car_option wealth_group)

gen co2_ann=prob_subs*2307.5*et_co2*average_use/1000
gen tax_total=prob_subs*2307.5*car_tax
gen tax_pv=prob_subs*2307.5*pv_tax
gen paid_subsidy=prob_subs*2307.5*subsidy if drivetype==3
gen emb_co2_ann=0
replace emb_co2_ann=prob_subs*2307.5*average_use/100*et_verbrauch*0.523 if drivetype==1
replace emb_co2_ann=prob_subs*2307.5*average_use/100*et_verbrauch*0.445 if drivetype==2
replace emb_co2_ann=prob_subs*2307.5*average_use/100*el_verbrauch*0.139 if drivetype==3
replace emb_co2_ann=prob_subs*2307.5*average_use/100*el_verbrauch*0.139 + prob_subs*2307.5*average_use/100*et_verbrauch*0.523 if drivetype==4


gen petrol_usage=prob_subs*2307.5*average_use/100*et_verbrauch
gen petrol_tax=petrol_usage*0.7312
replace petrol_tax=petrol_usage*0.7587 if drivetype==2
gen veh_tax=prob_subs*2307.5*0.04*prhnewprice
replace veh_tax=0 if drivetype==3


gcollapse (sum) prob_subs co2_ann emb_co2_ann petrol_tax tax_total tax_pv paid_subsidy veh_tax (firstnm) subsidy, by(drivetype wealth_group)

frame change cons`count'

gen CS_subs_change=CS_subs-CS_init
gen be_subs_change=be_subs-be_init
gcollapse (firstnm) CS_subs_change CS_init wealth_group subsidy (sum) be_subs_change be_init , by(id)
gcollapse (firstnm) subsidy (sum) CS_subs_change CS_init be_subs_change be_init, by(wealth_group)

gen CS_subs_pchange=100*CS_subs_change/abs(CS_init)
gen be_subs_pchange=100*be_subs_change/abs(be_init)

replace CS_subs_change=2307.5/5768.5*CS_subs_change
replace be_subs_change=2307.5/5768.5*be_subs_change

frame change default	

drop price_alt price_w_alt2 price_w_alt3 price_w_alt4 prob_subs CS_subs subsidy be_subs

local count=`count'+1

timer off 1
timer list 1

}


frame change frame1

forvalues ii=2/34 {
	
frameappend frame`ii', drop
}

frame change cons1

forvalues ii=2/34 {
	
frameappend cons`ii', drop
}



frame change frame1

gen prob_EV=prob_subs if drivetype==3
gen prob_HEV=prob_subs if drivetype==4
gen prob_gas=prob_subs if drivetype==1
gen prob_diesel=prob_subs if drivetype==2

gcollapse (sum) co2_ann emb_co2_ann petrol_tax tax_total tax_pv paid_subsidy veh_tax (firstnm) prob_EV prob_HEV prob_gas prob_diesel, by(subsidy wealth_group)

gen TOT_EINK=59514
replace TOT_EINK=87082 if wealth_group==2
replace TOT_EINK=115592 if wealth_group==3
replace TOT_EINK=195182 if wealth_group==4


frlink 1:1 wealth_group subsidy, frame(cons1) 

frget CS_subs_pchange CS_subs_change be_subs_change be_subs_pchange, from(cons1)

frame drop cons1

frame rename frame1 subs_sim

save subs_sim_new.dta, replace


replace co2_ann=co2_ann/1000
replace emb_co2_ann=emb_co2_ann/1000

** Create a graph of changes in probabilities:
frame put prob* subsidy wealth_group, into(prob_subs)
frame change prob_subs

collapse (mean) prob* , by(subsidy)

*Calculate probability changes:
frame create init_prob
frame change init_prob

import excel init_prob.xlsx, first clear

keep drivetype prob
gen obs=1

reshape wide prob, i(obs) j(drivetype)

rename prob1 prob_init_gas
rename prob2 prob_init_diesel
rename prob3 prob_init_EV
rename prob4 prob_init_HEV

frame change prob_subs
gen obs=1
frlink m:1 obs, frame(init_prob)

frget _all, from(init_prob)


gen prob_EV_d=100*(prob_EV-prob_init_EV)
gen prob_HEV_d=100*(prob_HEV-prob_init_HEV)
gen prob_gas_d=100*(prob_gas-prob_init_gas)
gen prob_diesel_d=100*(prob_diesel-prob_init_diesel)

*Graph the change in probabilities
twoway (line prob_EV_d subsidy, lpattern(solid) lwidth(medthick) lcolor(blue) ) (line prob_HEV_d subsidy, lpattern(dash) lwidth(medthick) lcolor(red)) (line prob_gas_d subsidy, lpattern(dot) lwidth(medthick) lcolor(green)) (line prob_diesel_d subsidy, lpattern(dash_dot) lwidth(medthick) lcolor(orange)), ytitle(Changes in adoption probabilities (p.p.)) xtitle(EV subsidy (CHF)) legend(order(1 "EV" 2 "Hybrid" 3 "Gasoline" 4 "Diesel")) clegend(span) graphregion(fcolor(white) lcolor(none)) plotregion(fcolor(white))

graph export "$root/Results/additional_results/Delta_probs_subs.png", replace

frame change subs_sim
frame drop prob_subs
frame drop init_prob

*Calculate the necessary steps to be able to create the graphs:
preserve
keep if subsidy==0
keep wealth_group co2_ann emb_co2_ann 
rename co2_ann co2_init
rename emb_co2_ann emb_co2_init
tempfile co2_init
save `co2_init'
restore
merge m:1 wealth_group using `co2_init'
drop _merge

gen d_co2_pipe=co2_ann-co2_init
gen d_co2_emb=emb_co2_ann-emb_co2_init
gen d_co2_tot=d_co2_pipe+d_co2_emb


gen co2_m=175*abs(d_co2_tot)

preserve
keep if subsidy==0
keep wealth_group tax_total
rename tax_total tax_init
tempfile tax_init
save `tax_init'
restore
merge m:1 wealth_group using `tax_init'
drop _merge

gen d_tax=tax_total-tax_init

preserve
keep if subsidy==0
keep wealth_group tax_pv
rename tax_pv tax_pv_init
tempfile tax_pv_init
save `tax_pv_init'
restore
merge m:1 wealth_group using `tax_pv_init'
drop _merge

gen d_tax_pv=tax_pv-tax_pv_init

preserve
keep if subsidy==0
keep wealth_group petrol_tax
rename petrol_tax petrol_tax_init
tempfile petrol_init
save `petrol_init'
restore
merge m:1 wealth_group using `petrol_init'
drop _merge

gen d_petrol_tax=petrol_tax-petrol_tax_init


preserve
keep if subsidy==0
keep wealth_group veh_tax
rename veh_tax veh_tax_init
tempfile veh_init
save `veh_init'
restore
merge m:1 wealth_group using `veh_init'
drop _merge

gen d_veh_tax=veh_tax-veh_tax_init



local r=0.06
local cap=1/`r' - (1/(`r'*(1+`r')^15))
di `cap'

gen d_pub=(-paid_subsidy+d_veh_tax+d_tax_pv+(d_petrol_tax)*`cap')/1000

gen d_wel=(CS_subs_change*1000-(be_subs_change*1000)-paid_subsidy+d_veh_tax+d_tax_pv+(abs(co2_m)+d_petrol_tax)*`cap')/1000
replace d_wel=0 if missing(d_wel)

*Generate the graph of emision reductions
sort subsidy wealth_group

twoway (line d_co2_tot subsidy if wealth_group==1, lpattern(solid) lwidth(medthick) lcolor(blue)) (line d_co2_tot subsidy if wealth_group==2, lpattern(dash) lwidth(medthick) lcolor(red)) (line d_co2_tot subsidy if wealth_group==3, lpattern(dot) lwidth(medthick) lcolor(green)) (line d_co2_tot subsidy if wealth_group==4, lpattern(dash_dot) lwidth(medthick) lcolor(orange)), ytitle(Change Emissions (tCO2)) xtitle(EV subsidy (CHF)) legend(order(1 "1st wealth quartile" 2 "2nd wealth quartile" 3 "3rd wealth quartile" 4 "4th wealth quartile")) clegend(span) graphregion(fcolor(white) lcolor(none)) plotregion(fcolor(white))
graph export "$root/Results/tables_graphs/appendix/EV_subs_co2_d.png", replace

*Create a graph of *experienced CS change: 
gen CSexp_subs_change=CS_subs_change-be_subs_change
twoway (line CSexp_subs_change subsidy if wealth_group==1, lpattern(solid) lwidth(medthick) lcolor(blue)) (line CSexp_subs_change subsidy if wealth_group==2, lpattern(dash) lwidth(medthick) lcolor(red)) (line CSexp_subs_change subsidy if wealth_group==3, lpattern(dot) lwidth(medthick) lcolor(green)) (line CSexp_subs_change subsidy if wealth_group==4, lpattern(dash_dot) lwidth(medthick) lcolor(orange)), ytitle(Consumer surplus (experienced) change (kCHF)) xtitle(EV subsidy (CHF)) legend(order(1 "1st wealth quartile" 2 "2nd wealth quartile" 3 "3rd wealth quartile" 4 "4th wealth quartile")) clegend(span) graphregion(fcolor(white) lcolor(none)) plotregion(fcolor(white))

graph export "$root/Results/tables_graphs/appendix/EV_subs-CS_change.png", replace

*Create a graph of change in public income
twoway (line d_pub subsidy if wealth_group==1, lpattern(solid) lwidth(medthick) lcolor(blue)) (line d_pub subsidy if wealth_group==2, lpattern(dash) lwidth(medthick) lcolor(red)) (line d_pub subsidy if wealth_group==3, lpattern(dot) lwidth(medthick) lcolor(green)) (line d_pub subsidy if wealth_group==4, lpattern(dash_dot) lwidth(medthick) lcolor(orange)), ytitle(Public finances change (kCHF)) xtitle(EV subsidy (CHF)) legend(order(1 "1st wealth quartile" 2 "2nd wealth quartile" 3 "3rd wealth quartile" 4 "4th wealth quartile")) clegend(span) graphregion(fcolor(white) lcolor(none)) plotregion(fcolor(white))

graph export "$root/Results/tables_graphs/appendix/EV_subs-pub_fin.png", replace

*Create a graph of the total welfare change: 

sort subsidy wealth_group
twoway (line d_wel subsidy if wealth_group==1, lpattern(solid) lwidth(medthick) lcolor(blue)) (line d_wel subsidy if wealth_group==2, lpattern(dash) lwidth(medthick) lcolor(red)) (line d_wel subsidy if wealth_group==3, lpattern(dot) lwidth(medthick) lcolor(green)) (line d_wel subsidy if wealth_group==4, lpattern(dash_dot) lwidth(medthick) lcolor(orange)), ytitle(Welfare change (kCHF)) xtitle(EV subsidy (CHF)) legend(order(1 "1st wealth quartile" 2 "2nd welth quartile" 3 "3rd wealth quartile" 4 "4th wealth quartile")) clegend(span) graphregion(fcolor(white) lcolor(none)) plotregion(fcolor(white))

graph export "$root/Results/tables_graphs/appendix/EV_subs_welfare.png", replace

frame change default
frame drop subs_sim


****
** Compute the third counterfactual - a vehicle registration tax related measure: 
****

frame change default

local r=0.06

gen pv_tax_alt=pv_tax

gen F_K=0.2
gen G_K=0.4

*for category F
replace pv_tax_alt=car_tax*(1+F_K)+(car_tax*(1+F_K))/((1+`r')^(1))+(car_tax*(1+F_K))/((1+`r')^(2))+(car_tax*(1+F_K))/((1+`r')^(3))+(car_tax)/((1+`r')^(4))+(car_tax)/((1+`r')^(5))+(car_tax)/((1+`r')^(6))+(car_tax)/((1+`r')^(7)) +(car_tax)/((1+`r')^(8))+(car_tax)/((1+`r')^(9))+(car_tax)/((1+`r')^(10))+(car_tax)/((1+`r')^(11))+(car_tax)/((1+`r')^(12))+(car_tax)/((1+`r')^(13))+(car_tax)/((1+`r')^(14)) if drivetype!=3 & env_cat==6
*for category G
replace pv_tax_alt=car_tax*(1+G_K)+(car_tax*(1+G_K))/((1+`r')^(1))+(car_tax*(1+G_K))/((1+`r')^(2))+(car_tax*(1+G_K))/((1+`r')^(3))+(car_tax)/((1+`r')^(4))+(car_tax)/((1+`r')^(5))+(car_tax)/((1+`r')^(6))+(car_tax)/((1+`r')^(7)) +(car_tax)/((1+`r')^(8))+(car_tax)/((1+`r')^(9))+(car_tax)/((1+`r')^(10))+(car_tax)/((1+`r')^(11))+(car_tax)/((1+`r')^(12))+(car_tax)/((1+`r')^(13))+(car_tax)/((1+`r')^(14)) if drivetype!=3 & env_cat==7


gen car_tax_alt=car_tax
replace car_tax_alt=car_tax*1.4 if drivetype!=3 & env_cat==7
replace car_tax_alt=car_tax*1.2 if drivetype!=3 & env_cat==6


gen alt_drive_cost=drive_cost_ann+drive_cost_ann*((1-((1+`r')^(-14)))/`r')

gen alt_var_cost=(alt_drive_cost+pv_tax_alt)/1000

replace alt_var_cost=0 if missing(alt_var_cost)

drop alt_drive_cost  F_K G_K

sort id car_option


mata
X=. 
st_view(X,.,"price height weight alt_var_cost typkw ev hyrid diesel car_size size_hh2 size_hh3 size_hh4 size_hh5 age_kw2 age_kw3 ev_urban2 ev_urban3 ev_dist_ev ev_range5 ev_year2 ev_year3 ev_pv ev_ho ev_wealth2 ev_wealth3 ev_wealth4 brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE2 car_cat_FE3 car_cat_FE4 car_cat_FE5 car_cat_FE6 car_cat_FE7 car_cat_FE8 env_friendly resid price_wealth2 price_wealth3 price_wealth4") 

l_fit=exp(X*beta)

var_cost=. ;
st_view(var_cost,.,"alt_var_cost") ;


/*
for (i=1; i<=ndraws; i++) {
st_addvar("double", "fit_" + strofreal(i))
st_store(.,"fit_" + strofreal(i),l_fit[,i])
}
*/
id=st_data(.,"id")
V=panelsetup(id, 1)

    ncol   = ndraws
    bigtot = J(rows(id), ncol, .)
    for (i = 1; i <= rows(V); i++) {
        X1 = panelsubmatrix(l_fit, i, V)
        bigtot[|V[i,1], 1 \ V[i, 2], ncol|] = J(rows(X1), 1, colsum(X1))
}

mata drop X1 ncol V

_editmissing(l_fit,0)

prob_all=l_fit:/bigtot
mata drop l_fit X 

prob=mean(prob_all')

be=prob_all :* (J(1, cols(gamma), var_cost)-(gamma:*var_cost)) ;

be_m=mean(be') ;

mata drop prob_all be

st_addvar("double", "prob_count3")
st_store(.,"prob_count3",prob')

st_addvar("double", "be_count3") ;
st_store(.,"be_count3",be_m') ;



CS=mean((((-1):/b_price):*log(bigtot))')


mata drop bigtot 

st_addvar("double", "CS_count3")
st_store(.,"CS_count3",CS')


mata drop prob CS be_m var_cost
end


*First calculate the predicted probabilites for each category:
frame put prob_count3 et_co2 average_use et_verbrauch el_verbrauch car_tax_alt pv_tax_alt drivetype car_option prhnewprice, into(alt3_prob)
frame change alt3_prob

gcollapse (mean) prob_count3 et_co2 average_use et_verbrauch el_verbrauch car_tax_alt pv_tax_alt prhnewprice (firstnm) drivetype , by(car_option)


*We use the number of 9230 cars as the average registrations per year

gen co2_ann=prob_count3*9230*et_co2*average_use/1000
gen emb_co2_ann=0
replace emb_co2_ann=prob*9230*average_use/100*et_verbrauch*0.523 if drivetype==1
replace emb_co2_ann=prob*9230*average_use/100*et_verbrauch*0.445 if drivetype==2
replace emb_co2_ann=prob*9230*average_use/100*el_verbrauch*0.139 if drivetype==3
replace emb_co2_ann=prob*9230*average_use/100*el_verbrauch*0.139 + prob*9230*average_use/100*et_verbrauch*0.523 if drivetype==4


gen petrol_usage=prob*9230*average_use/100*et_verbrauch
gen petrol_tax=petrol_usage*0.7312
replace petrol_tax=petrol_usage*0.7587 if drivetype==2

gen tax_total=prob_count3*9230*car_tax_alt
gen tax_pv=prob_count3*9230*pv_tax_alt
gen veh_tax=prob_count3*9230*prhnewprice*0.04
replace veh_tax=0 if drivetype==3

gcollapse (sum) prob_count3 co2_ann emb_co2_ann petrol_tax  tax_total tax_pv veh_tax, by(drivetype)

export excel "alt3_prob.xlsx", firstrow(variables) replace

frame change default
frame drop alt3_prob

*Also calculate the predicted probabilites for each category by wealth quartile

frame put prob_count3 et_co2 average_use et_verbrauch el_verbrauch car_tax_alt pv_tax_alt drivetype car_option wealth_group prhnewprice, into(alt3_byinc)

frame change alt3_byinc

gcollapse (mean) prob_count3 et_co2 average_use et_verbrauch el_verbrauch car_tax_alt pv_tax_alt prhnewprice (firstnm) drivetype, by(car_option wealth_group)

gen co2_ann=prob_count3*2307.5*et_co2*average_use/1000
gen tax_total=prob_count3*2307.5*car_tax_alt
gen tax_pv=prob_count3*2307.5*pv_tax_alt
gen emb_co2_ann=0
replace emb_co2_ann=prob*2307.5*average_use/100*et_verbrauch*0.523 if drivetype==1
replace emb_co2_ann=prob*2307.5*average_use/100*et_verbrauch*0.445 if drivetype==2
replace emb_co2_ann=prob*2307.5*average_use/100*el_verbrauch*0.139 if drivetype==3
replace emb_co2_ann=prob*2307.5*average_use/100*el_verbrauch*0.139 + prob*2307.5*average_use/100*et_verbrauch*0.523 if drivetype==4


gen petrol_usage=prob*2307.5*average_use/100*et_verbrauch
gen petrol_tax=petrol_usage*0.7312
replace petrol_tax=petrol_usage*0.7587 if drivetype==2
gen veh_tax=prob_count3*2307.5*prhnewprice*0.04
replace veh_tax=0 if drivetype==3


collapse (sum) prob_count3 co2_ann emb_co2_ann petrol_tax  tax_total tax_pv veh_tax, by(drivetype wealth_group)

export excel "alt3_prob_byinc.xlsx", firstrow(variables) replace

frame change default
frame drop alt3_byinc


**Calculate the change in consumer welfare 


*Now we can calculate consumer surplus based on the adequate formula of logsum
*check which income was used in formula for log(price)/inc

gen CS_change3=CS_count3-CS_init
gen be_change3=be_count3-be_init


frame put CS_change3 CS_init wealth_group id be_change3 be_init, into(cs3)
frame change cs3
rename CS_change3 CS_change
rename be_change3 be_change
collapse (firstnm) CS_change CS_init wealth_group (sum) be_change be_init, by(id)
collapse (sum) CS_change CS_init be_init be_change , by(wealth_group)
gen CS_perc_change=CS_change/abs(CS_init)
gen be_perc_change=be_change/abs(be_init)

replace CS_change=2307.5/5768.5*CS_change
replace be_change=2307.5/5768.5*be_change


export excel "CS3.xlsx", firstrow(variables) replace

frame change default
frame drop cs3



**Add an additional program here where a varying tax break / feebate is calculated - First need to discuss how this should look!

*At level zero, the results represent the initial probability analysis: 
frame put prob et_co2 average_use et_verbrauch el_verbrauch car_tax pv_tax drivetype car_option wealth_group prhnewprice, into(init)

frame change init

rename prob prob_taxsim
gen tax_sim=0

gcollapse (mean) prob_taxsim et_co2 average_use et_verbrauch el_verbrauch car_tax pv_tax prhnewprice (firstnm) drivetype tax_sim, by(car_option wealth_group)

gen co2_ann=prob_taxsim*2307.5*et_co2*average_use/1000
gen tax_total=prob_taxsim*2307.5*car_tax
gen tax_pv=prob_taxsim*2307.5*pv_tax
gen emb_co2_ann=0
replace emb_co2_ann=prob_taxsim*2307.5*average_use/100*et_verbrauch*0.523 if drivetype==1
replace emb_co2_ann=prob_taxsim*2307.5*average_use/100*et_verbrauch*0.445 if drivetype==2
replace emb_co2_ann=prob_taxsim*2307.5*average_use/100*el_verbrauch*0.139 if drivetype==3
replace emb_co2_ann=prob_taxsim*2307.5*average_use/100*el_verbrauch*0.139 + prob_taxsim*2307.5*average_use/100*et_verbrauch*0.523 if drivetype==4


gen petrol_usage=prob*2307.5*average_use/100*et_verbrauch
gen petrol_tax=petrol_usage*0.7312
replace petrol_tax=petrol_usage*0.7587 if drivetype==2
gen veh_tax=prob_taxsim*2307.5*prhnewprice*0.04
replace veh_tax=0 if drivetype==3



gcollapse (sum) prob_taxsim co2_ann emb_co2_ann tax_total tax_pv petrol_tax veh_tax (firstnm) tax_sim, by(drivetype wealth_group)


frame put _all, into(frame1)

frame change default
frame drop init



*At level zero the CS change is also zero
frame create cons1
frame change cons1

set obs 1
gen CS_tax_pchange=0
gen be_tax_pchange=0
gen CS_tax_change=0
gen be_tax_change=0

gen tax_sim=0


frame change default
local count=2

drop car_tax_alt pv_tax_alt

forvalues x=1/14 {
gen car_tax_alt=car_tax
replace car_tax_alt=car_tax/0.4 if drivetype==3
replace car_tax_alt=car_tax/0.6 if drivetype!=3 & env_cat==1
replace car_tax_alt=car_tax/0.8 if drivetype!=3 & env_cat==2	
	
	gen tax_sim=`x'
	
*Add the tax feebates:
	merge m:1 tax_sim using "$root/Data/Produced/tax_params.dta", 
	drop if _merge==2
	drop _merge
	
	
	
	
local r=0.06

gen pv_tax_alt=pv_tax

*for EV
replace pv_tax_alt=car_tax_alt*(1+EV_K)+(car_tax_alt*(1+EV_K))/((1+`r')^(1))+(car_tax_alt*(1+EV_K))/((1+`r')^(2))+(car_tax_alt*(1+EV_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9))+(car_tax_alt)/((1+`r')^(10))+(car_tax_alt)/((1+`r')^(11))+(car_tax_alt)/((1+`r')^(12))+(car_tax_alt)/((1+`r')^(13))+(car_tax_alt)/((1+`r')^(14)) if drivetype==3 

*for category A
replace pv_tax_alt=car_tax_alt*(1+A_K)+(car_tax_alt*(1+A_K))/((1+`r')^(1))+(car_tax_alt*(1+A_K))/((1+`r')^(2))+(car_tax_alt*(1+A_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9))+(car_tax_alt)/((1+`r')^(10))+(car_tax_alt)/((1+`r')^(11))+(car_tax_alt)/((1+`r')^(12))+(car_tax_alt)/((1+`r')^(13))+(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==1

*for category B
replace pv_tax_alt=car_tax_alt*(1+B_K)+(car_tax_alt*(1+B_K))/((1+`r')^(1))+(car_tax_alt*(1+B_K))/((1+`r')^(2))+(car_tax_alt*(1+B_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9))+(car_tax_alt)/((1+`r')^(10))+(car_tax_alt)/((1+`r')^(11))+(car_tax_alt)/((1+`r')^(12))+(car_tax_alt)/((1+`r')^(13))+(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==2

*for category E
replace pv_tax_alt=car_tax_alt*(1+E_K)+(car_tax_alt*(1+E_K))/((1+`r')^(1))+(car_tax_alt*(1+E_K))/((1+`r')^(2))+(car_tax_alt*(1+E_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9))+(car_tax_alt)/((1+`r')^(10))+(car_tax_alt)/((1+`r')^(11))+(car_tax_alt)/((1+`r')^(12))+(car_tax_alt)/((1+`r')^(13))+(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==5
*for category F
replace pv_tax_alt=car_tax_alt*(1+F_K)+(car_tax_alt*(1+F_K))/((1+`r')^(1))+(car_tax_alt*(1+F_K))/((1+`r')^(2))+(car_tax_alt*(1+F_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9))+(car_tax_alt)/((1+`r')^(10))+(car_tax_alt)/((1+`r')^(11))+(car_tax_alt)/((1+`r')^(12))+(car_tax_alt)/((1+`r')^(13))+(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==6
*for category G
replace pv_tax_alt=car_tax_alt*(1+G_K)+(car_tax_alt*(1+G_K))/((1+`r')^(1))+(car_tax_alt*(1+G_K))/((1+`r')^(2))+(car_tax_alt*(1+G_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9))+(car_tax_alt)/((1+`r')^(10))+(car_tax_alt)/((1+`r')^(11))+(car_tax_alt)/((1+`r')^(12))+(car_tax_alt)/((1+`r')^(13))+(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==7

replace car_tax_alt=car_tax_alt*(1+G_K) if drivetype!=3 & env_cat==7
replace car_tax_alt=car_tax_alt*(1+F_K) if drivetype!=3 & env_cat==6
replace car_tax_alt=car_tax_alt*(1+E_K) if drivetype!=3 & env_cat==5
replace car_tax_alt=car_tax_alt*(1+B_K) if drivetype!=3 & env_cat==2
replace car_tax_alt=car_tax_alt*(1+A_K) if drivetype!=3 & env_cat==1
replace car_tax_alt=car_tax_alt*(1+EV_K) if drivetype==3 


gen alt_drive_cost=drive_cost_ann+drive_cost_ann*((1-((1+`r')^(-14)))/`r')

gen alt_var_cost4=(alt_drive_cost+pv_tax_alt)/1000

replace alt_var_cost4=0 if missing(alt_var_cost4)

drop alt_drive_cost 
	
** Add here probability calculation:

sort id car_option


mata {
X=. ;
st_view(X,.,"price height weight alt_var_cost4 typkw ev hyrid diesel car_size size_hh2 size_hh3 size_hh4 size_hh5 age_kw2 age_kw3 ev_urban2 ev_urban3 ev_dist_ev ev_range5 ev_year2 ev_year3 ev_pv ev_ho ev_wealth2 ev_wealth3 ev_wealth4 brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE2 car_cat_FE3 car_cat_FE4 car_cat_FE5 car_cat_FE6 car_cat_FE7 car_cat_FE8 env_friendly resid price_wealth2 price_wealth3 price_wealth4") ;


l_fit=exp(X*beta) ;


var_cost=. ;
st_view(var_cost,.,"alt_var_cost4") ;

/*
for (i=1; i<=ndraws; i++) {
st_addvar("double", "fit_" + strofreal(i))
st_store(.,"fit_" + strofreal(i),l_fit[,i])
}
*/
id=st_data(.,"id") ;
V=panelsetup(id, 1) ;

    ncol   = ndraws ;
    bigtot = J(rows(id), ncol, .) ;
    for (i = 1; i <= rows(V); i++) { ;
        X1 = panelsubmatrix(l_fit, i, V) ;
        bigtot[|V[i,1], 1 \ V[i, 2], ncol|] = J(rows(X1), 1, colsum(X1)) ;
} ;
/* instead of dropping recode the memory intensive part
mata drop X1 ncol V ;
*/
X1=. ;
V=. ;
_editmissing(l_fit,0) ;

prob_all=l_fit:/bigtot ;


/*
mata drop l_fit X  ;
*/
l_fit=.
X=.
prob=mean(prob_all') ;

be=prob_all :* (J(1, cols(gamma), var_cost)-(gamma:*var_cost)) ;

be_m=mean(be') ;

/*
mata drop prob_all ;
*/
prob_all=. ;
be=. ;
var_cost=.;
st_addvar("double", "prob_taxsim") ;
st_store(.,"prob_taxsim",prob') ;

st_addvar("double", "be_taxsim") ;
st_store(.,"be_taxsim",be_m') ;


CS=mean((((-1):/b_price):*log(bigtot))') ;

/*
mata drop bigtot ;
*/
bigtot=.
st_addvar("double", "CS_taxsim") ;
st_store(.,"CS_taxsim",CS') ;

/*
mata drop prob CS ;
*/
}


*Create a frame to calculate the consumer welfare change:
frame put CS_taxsim CS_init be_taxsim be_init id wealth_group tax_sim, into(cons`count')


frame put tax_sim prob_taxsim et_co2 average_use et_verbrauch el_verbrauch car_tax_alt pv_tax_alt drivetype car_option wealth_group prhnewprice, into(frame`count')

frame change frame`count'

gcollapse (mean) prob_taxsim et_co2 average_use et_verbrauch el_verbrauch car_tax_alt pv_tax_alt prhnewprice (firstnm) drivetype tax_sim, by(car_option wealth_group)

gen co2_ann=prob_taxsim*2307.5*et_co2*average_use/1000
gen tax_total=prob_taxsim*2307.5*car_tax_alt
gen tax_pv=prob_taxsim*2307.5*pv_tax_alt
gen emb_co2_ann=0
replace emb_co2_ann=prob_taxsim*2307.5*average_use/100*et_verbrauch*0.523 if drivetype==1
replace emb_co2_ann=prob_taxsim*2307.5*average_use/100*et_verbrauch*0.445 if drivetype==2
replace emb_co2_ann=prob_taxsim*2307.5*average_use/100*el_verbrauch*0.139 if drivetype==3
replace emb_co2_ann=prob_taxsim*2307.5*average_use/100*el_verbrauch*0.139 + prob_taxsim*2307.5*average_use/100*et_verbrauch*0.523 if drivetype==4


gen petrol_usage=prob_taxsim*2307.5*average_use/100*et_verbrauch
gen petrol_tax=petrol_usage*0.7312
replace petrol_tax=petrol_usage*0.7587 if drivetype==2
gen veh_tax=prob_taxsim*2307.5*prhnewprice*0.04
replace veh_tax=0 if drivetype==3



gcollapse (sum) prob_taxsim co2_ann emb_co2_ann petrol_tax tax_total tax_pv veh_tax (firstnm) tax_sim, by(drivetype wealth_group)

frame change cons`count'

gcollapse (firstnm) CS_taxsim CS_init wealth_group tax_sim (sum) be_init be_taxsim, by(id)
gen CS_tax_change=CS_taxsim-CS_init
gen be_tax_change=be_taxsim-be_init
gcollapse (firstnm) tax_sim (sum) CS_tax_change CS_init be_tax_change be_init, by(wealth_group)

replace CS_tax_change=2307.5/5768.5*CS_tax_change
replace be_tax_change=2307.5/5768.5*be_tax_change

gen CS_tax_pchange=100*CS_tax_change/abs(CS_init)
gen be_tax_pchange=100*be_tax_change/abs(be_init)

frame change default	

local count=`count'+1

drop prob_taxsim tax_sim CS_taxsim alt_var_cost4 car_tax_alt pv_tax_alt A_K B_K E_K F_K G_K EV_K be_taxsim

}



frame change frame1

forvalues ii=2/15 {
	
frameappend frame`ii', drop
}

frame change cons1

forvalues ii=2/15 {
	
frameappend cons`ii', drop
}


frame change frame1


gen prob_EV=prob_taxsim if drivetype==3
gen prob_HEV=prob_taxsim if drivetype==4
gen prob_gas=prob_taxsim if drivetype==1
gen prob_diesel=prob_taxsim if drivetype==2

collapse (sum) co2_ann emb_co2_ann petrol_tax tax_total tax_pv veh_tax (firstnm) prob_EV prob_HEV prob_gas prob_diesel , by(tax_sim wealth_group)

gen TOT_EINK=59514
replace TOT_EINK=87082 if wealth_group==2
replace TOT_EINK=115592 if wealth_group==3
replace TOT_EINK=195182 if wealth_group==4


gen tax_incid=tax_total/(2307.5*TOT_EINK)

frame rename frame1 tax_sim


frlink 1:1 wealth_group tax_sim, frame(cons1) 

frget CS_tax_pchange CS_tax_change be_tax_change be_tax_pchange, from(cons1)

frame drop cons1

save tax_sim_new.dta, replace



replace co2_ann=co2_ann/1000
replace emb_co2_ann=emb_co2_ann/1000


** Create a graph of changes in probabilities:
frame put prob* tax_sim wealth_group, into(prob_tax)
frame change prob_tax

collapse (mean) prob* , by(tax_sim)

*Calculate probability changes:
frame create init_prob
frame change init_prob

import excel init_prob.xlsx, first clear

keep drivetype prob
gen obs=1

reshape wide prob, i(obs) j(drivetype)

rename prob1 prob_init_gas
rename prob2 prob_init_diesel
rename prob3 prob_init_EV
rename prob4 prob_init_HEV

frame change prob_tax
gen obs=1
frlink m:1 obs, frame(init_prob)

frget _all, from(init_prob)

*Calculate probability changes:

gen prob_EV_d=100*(prob_EV-prob_init_EV)
gen prob_HEV_d=100*(prob_HEV-prob_init_HEV)
gen prob_gas_d=100*(prob_gas-prob_init_gas)
gen prob_diesel_d=100*(prob_diesel-prob_init_diesel)

*Graph the change in probabilities
twoway (line prob_EV_d tax_sim if tax_sim>0, lpattern(solid) lwidth(medthick) lcolor(blue) ) (line prob_HEV_d tax_sim if tax_sim>0, lpattern(dash_dot) lwidth(medthick) lcolor(orange)) (line prob_gas_d tax_sim if tax_sim>0, lpattern(dash) lwidth(medthick) lcolor(red)) (line prob_diesel_d tax_sim if tax_sim>0, lpattern(dot) lwidth(medthick) lcolor(green)) , ytitle(Changes in adoption probabilities (p.p.)) xtitle(Vehicle registration tax simulation) legend(order(1 "EV" 2 "Hybrid" 3 "Gasoline" 4 "Diesel")) clegend(span) graphregion(fcolor(white) lcolor(none)) plotregion(fcolor(white))

graph export "$root/Results/additional_results/Delta_probs_tax.png" , replace

frame change tax_sim
frame drop prob_tax
frame drop init_prob

*Calculate the necessary results to create the graph:
preserve
keep if tax_sim==0
keep wealth_group co2_ann emb_co2_ann
rename co2_ann co2_init
rename emb_co2_ann emb_co2_ann_init
tempfile co2_init
save `co2_init'
restore
merge m:1 wealth_group using `co2_init'
drop _merge

gen d_co2=co2_ann-co2_init
gen d_emb_co2=emb_co2_ann-emb_co2_ann_init
gen d_co2_tot=d_co2+d_emb_co2

*Calculate total welfare change: 
gen co2_m=175*abs(d_co2_tot)

preserve
keep if tax_sim==0
keep wealth_group tax_total
rename tax_total tax_init
tempfile tax_init
save `tax_init'
restore
merge m:1 wealth_group using `tax_init'
drop _merge

gen d_tax=tax_total-tax_init

preserve
keep if tax_sim==0
keep wealth_group tax_pv
rename tax_pv tax_pv_init
tempfile tax_pv_init
save `tax_pv_init'
restore
merge m:1 wealth_group using `tax_pv_init'
drop _merge

gen d_tax_pv=tax_pv-tax_pv_init

preserve
keep if tax_sim==0
keep wealth_group petrol_tax
rename petrol_tax petrol_tax_init
tempfile petrol_init
save `petrol_init'
restore
merge m:1 wealth_group using `petrol_init'
drop _merge

gen d_petrol_tax=petrol_tax-petrol_tax_init


preserve
keep if tax_sim==0
keep wealth_group veh_tax
rename veh_tax veh_tax_init
tempfile veh_init
save `veh_init'
restore
merge m:1 wealth_group using `veh_init'
drop _merge

gen d_veh_tax=veh_tax-veh_tax_init

local r=0.06
local cap=1/`r' - (1/(`r'*(1+`r')^15))
di `cap'

gen d_pub=(d_veh_tax + d_tax_pv + `cap'*(d_petrol_tax))/1000

gen d_wel=(CS_tax_change*1000-(be_tax_change*1000)+d_veh_tax+d_tax_pv + (abs(co2_m)+d_petrol_tax)*`cap')/1000
replace d_wel=0 if missing(d_wel)


*Emission reduction graph:
sort tax_sim wealth_group

twoway (line d_co2_tot tax_sim if wealth_group==1 & tax_sim>0, lpattern(solid) lwidth(medthick) lcolor(blue)) (line d_co2_tot tax_sim if wealth_group==2 & tax_sim>0, lpattern(dash) lwidth(medthick) lcolor(red)) (line d_co2_tot tax_sim if wealth_group==3 & tax_sim>0, lpattern(dot) lwidth(medthick) lcolor(red)) (line d_co2_tot tax_sim if wealth_group==4 & tax_sim>0, lpattern(dash_dot) lwidth(medthick) lcolor(orange)), ytitle(Change Emissions (tCO2)) xtitle(Vehicle registration tax simulation) legend(order(1 "1st wealth quartile" 2 "2nd wealth quartile" 3 "3rd wealth quartile" 4 "4th wealth quartile")) clegend(span) graphregion(fcolor(white) lcolor(none)) plotregion(fcolor(white))
graph export "$root/Results/tables_graphs/appendix/tax_sim_co2_d.png", replace


*Generate a graph of change in public finance: 
twoway (line d_pub tax_sim if wealth_group==1 & tax_sim>0, lpattern(solid) lwidth(medthick) lcolor(blue)) (line d_pub tax_sim if wealth_group==2 & tax_sim>0, lpattern(dash) lwidth(medthick) lcolor(red)) (line d_pub tax_sim if wealth_group==3 & tax_sim>0, lpattern(dot) lwidth(medthick) lcolor(green)) (line d_pub tax_sim if wealth_group==4 & tax_sim>0, lpattern(dash_dot) lwidth(medthick) lcolor(orange)), ytitle(Change public finance (kCHF NPV)) xtitle(Vehicle registration tax simulation) legend(order(1 "1st wealth quartile" 2 "2nd wealth quartile" 3 "3rd wealth quartile" 4 "4th wealth quartile")) clegend(span) graphregion(fcolor(white) lcolor(none)) plotregion(fcolor(white))

graph export "$root/Results/tables_graphs/appendix/tax_sim-d_pub.png", replace

*Generate a graph of Consumer surplus change
gen CSexp_tax_change=CS_tax_change-be_tax_change
twoway (line CSexp_tax_change tax_sim if wealth_group==1 & tax_sim>0, lpattern(solid) lwidth(medthick) lcolor(blue)) (line CSexp_tax_change tax_sim if wealth_group==2 & tax_sim>0, lpattern(dash) lwidth(medthick) lcolor(red)) (line CSexp_tax_change tax_sim if wealth_group==3 & tax_sim>0, lpattern(dot) lwidth(medthick) lcolor(green)) (line CSexp_tax_change tax_sim if wealth_group==4 & tax_sim>0, lpattern(dash_dot) lwidth(medthick) lcolor(orange)), ytitle(Consumer surplus (experienced) change (kCHF)) xtitle(Vehicle registration tax simulation) legend(order(1 "1st wealth quartile" 2 "2nd wealth quartile" 3 "3rd wealth quartile" 4 "4th wealth quartile")) clegend(span) graphregion(fcolor(white) lcolor(none)) plotregion(fcolor(white))

graph export "$root/Results/tables_graphs/appendix/tax_sim-CS_d.png", replace
*Welfare change:
sort tax_sim wealth_group
twoway (line d_wel tax_sim if wealth_group==1 & tax_sim>0, lpattern(solid) lwidth(medthick) lcolor(blue)) (line d_wel tax_sim if wealth_group==2 & tax_sim>0, lpattern(dash) lwidth(medthick) lcolor(red)) (line d_wel tax_sim if wealth_group==3 & tax_sim>0, lpattern(dot) lwidth(medthick) lcolor(green)) (line d_wel tax_sim if wealth_group==4 & tax_sim>0, lpattern(dash_dot) lwidth(medthick) lcolor(orange)), ytitle(Welfare change (kCHF)) xtitle(Vehicle registration tax simulation) legend(order(1 "1st wealth quartile" 2 "2nd wealth quartile" 3 "3rd wealth quartile" 4 "4th wealth quartile")) clegend(span) graphregion(fcolor(white) lcolor(none)) plotregion(fcolor(white))

graph export "$root/Results/tables_graphs/appendix/tax_sim_welfare.png", replace


 
frame change default
frame drop tax_sim


